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A  finite  strain  theory  is  developed  for  anisotropic  single  crystals  undergoing  shock  loading. 
Inelastic  deformation  may  arise  from  dislocation  slip,  twinning,  or  fracture  and  crack  slid¬ 
ing.  Internal  energy  can  generally  depend  on  a  logarithmic  measure  of  finite  elastic  strain, 
entropy,  and  an  internal  variable  associated  with  defect  accumulation.  A  closed  form  ana¬ 
lytical  solution  is  derived  for  the  planar  shock  response  in  the  thermoelastic  regime,  at 
axial  stresses  up  to  the  Hugoniot  elastic  limit.  In  the  plastic  regime,  for  highly  symmetric 
orientations  and  rate  independent  shear  strength,  the  Rankine-Hugoniot  conditions  and 
constitutive  relations  can  be  reduced  to  a  set  of  algebraic  equations  that  can  be  solved 
for  the  material  response.  The  theory  is  applied  towards  planar  shock  loading  of  single 
crystals  of  sapphire,  diamond,  and  quartz.  Logarithmic  elasticity  is  demonstrated  to  be 
more  accurate  (i.e.,  require  fewer  higher-order  elastic  constants)  than  Lagrangian  or  Eule- 
rian  theories  for  sapphire,  diamond,  and  Z-cut  quartz.  Results  provide  new  insight  into  cri¬ 
teria  for  initiation  of  twinning,  slip,  and/or  fracture  in  these  materials  as  well  as  their 
strength  degradation  when  shocked  at  increasingly  higher  pressures  above  the  Hugoniot 
elastic  limit. 
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1.  Introduction 

The  response  of  solids  to  shock  compression  under  planar  impact  has  been  a  subject  of  intensive  study  over  the  past  half- 
century  (McQueen,  Marsh,  Taylor,  Fritz,  &  Carter,  1970),  including  numerous  advances  in  experiments,  theoretical/analytical 
methods,  and  numerical  techniques.  Regarding  modeling  of  related  phenomena,  much  effort  has  centered  on  development  of 
the  pressure-volume  equation-of-state  (EOS)  applicable  for  loading  regimes  or  materials  (e.g.,  very  high  pressures,  or  isotro¬ 
pic  fluids  and  ductile  solids)  wherein  shear  strength  and  anisotropy  are  of  little  or  no  concern.  However,  for  strong  solids- 
ceramics,  minerals,  and  some  metals  and  alloys,  for  example  -  significant  shear  strength  is  retained  under  impact  loading  at 
moderate  to  high  pressures.  This  strength  can  affect  the  global  response  of  the  material  in  loading  conditions  pertinent  to 
ballistic  penetration,  geologic  events,  explosions,  high-speed  vehicular  collisions,  etc.  Microstructure,  including  grain  or  con¬ 
stituent  orientation  and  presence  of  multiple  phases,  can  also  significantly  affect  the  shock  response  (Grady,  1984).  In 
shocked  single  crystals  of  high  purity  which  are  the  focus  of  the  present  work,  crystal  lattice  orientation  is  the  primary 
descriptor  of  initial  microstructure,  and  it  affects  anisotropic  themoelasticity  and  orientation-dependent  inelasticity  (e.g., 
slip,  twinning,  and  cleavage  fracture). 
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In  Section  2  of  this  paper,  a  new  finite  deformation,  anisotropic  thermoelastic  theory  is  developed  for  single  crystals,  and 
is  applied  to  study  the  shock  response  of  oriented  crystals  of  sapphire  (oc-A/203),  diamond  (C),  and  quartz  (a-Si02).  These 
materials  are  considered  because  (i)  they  exhibit  a  high  Hugoniot  Elastic  Limit  (HEL),  enabling  assessment  of  finite  strain 
effects  in  their  elastic  shock  response  and  (ii)  ample  data  on  higher-order  elastic  constants  (Graham,  1972;  Hankey  & 
Schuele,  1970;  Nielsen,  1986;  Thurston,  McSkimin,  &  Andreatch,  1966)  and  planar  shock  compression  experiments  (Graham 
&  Brooks,  1971;  Fowles,  1967;  Lang  &  Gupta,  2010)  exist,  enabling  model  development  and  validation. 

In  recent  prior  work  (Clayton,  2013),  a  finite  strain  theory  based  on  strain  measure  D  =1(1  -F~'F~t)  was  developed  and 
compared  to  usual  nonlinear  crystal  thermoelasticity  (Clayton,  2011a;  Thurston,  1974;  Wallace,  1972)  based  on  the  Green 
strain  E  =  \  (FTF  -  1),  where  F  is  the  deformation  gradient.  These  are  respectively  referred  to  as  “Eulerian”  and  “Lagrangian" 
theories  (Nielsen,  1986;  Thomsen,  1972),  though  both  strain  measures  are  referred  to  material  coordinates  and  are  thus 
admissible  in  thermodynamic  potentials  for  anisotropic  solids  [in  contrast  to  the  Almansi  strain  A  =  1(1  -  F~TF~'),  for  exam¬ 
ple,  used  in  isotropic  nonlinear  elasticity].  Eulerian  theory  demonstrated  advantages  with  regards  to  describing  the  compres¬ 
sion  and  shear  reponses  of  ideal  cubic  crystals  with  Cauchy  symmetry,  and  demonstrated  greater  intrinsic  stability  (Clayton 
&  Bliss,  2014)  than  Lagrangian  theory  for  these  conditions.  Analytical  solutions  to  the  planar  shock  problem  were  also  de¬ 
rived  using  both  models,  with  predictions  compared  for  sapphire,  diamond,  and  quartz  (Clayton,  2013).  Neither  model 
was  definitively  superior  for  describing  the  shock  response  of  these  real  anisotropic  single  crystals,  with  elastic  constants 
of  up  to  order  three  or  four  necessary  in  either  case.  Anisotropic  Lagrangian-type  nonlinear  elastic  models  of  crystals  have 
been  used  by  other  authors  in  numerical  simulations  of  wave  propagation  (Winey  &  Gupta,  2004)  and  spall  (Foulk  8;  Vogler, 
2010).  Recently,  the  aforementioned  Eulerian  theory  has  been  applied  towards  new  nonlinear  elastic  solutions  of  boundary 
value  problems  involving  discrete  lattice  defects  (Clayton,  in  press). 

In  the  novel  thermoelastic  theory  developed  in  the  current  work,  internal  energy  is  a  function  of  entropy  and  material 
logarithmic  strain  e  =  In  U,  where  U  is  the  right  stretch  in  the  polar  decomposition  F  =  RU  =  VR,  with  R  the  rotation.  Elastic 
theory  based  on  Hencky’s  strain  measure  in  V  has  been  used  to  accurately  model  isotropic  solids  at  moderate-to-large  strains 
(Anand,  1979),  but  this  Eulerian  theory  does  not  apply  for  anisotropic  solids.  Regarding  the  latter,  hyperelastic  theory  based 
on  e  has  been  considered  with  regards  to  derivation  of  higher-order  elastic  constants  in  cubic  crystals  (Dluzewski,  2000),  but 
such  theory  has  remained,  until  now,  untested  for  stress  states  involving  both  pressure  and  shear  such  as  uniaxial  strain 
shock  compression.  Success  of  the  logarithmic  pressure-volume  EOS,  to  which  e-based  theory  degenerates  under  hydrostatic 
loading,  has  been  demonstrated  (Poirier  &  Tarantola,  1998).  In  Section  2.1,  a  complete  thermodynamically  consistent  non¬ 
linear  thermoelasticity  theory  incorporating  e  is  derived  and  presented  for  the  first  time.  A  new  analytical  solution  to  the 
planar  shock  problem  for  solids  obeying  this  constitutive  theory  is  derived  in  Section  2.2.  Application  of  the  solution  to  sap¬ 
phire,  diamond,  and  quartz  follows  in  Section  2.3,  wherein  advantages  of  the  proposed  logarithmic  formulation  over  existing 
Lagrangian  and  Eulerian  thermoelasticity  models  become  clear. 

Recently  (Srinivasa,  2012),  a  promising  structure  for  anisotropic  hyperelasticity  has  been  developed  that  involves  decom¬ 
position  of  F  into  the  product  of  an  orthogonal  matrix  and  an  upper  triangular  matrix,  with  strain  energy  depending  on  the 
latter.  This  approach,  though  not  explored  in  the  current  paper,  demonstrates  certain  advantages  regarding  computational 
efficiency  and  physical  interpretation  relative  to  models  that  use  a  polar  decomposition  (e.g.,  logarithmic  theory). 

In  Section  3,  the  logarithmic  theory  is  extended  to  address  inelastic  deformation,  wherein  the  deformation  gradient  is 
split  into  thermoelastic  (FE)  and  inelastic  (Fp)  parts:  F  =  FEFP  (Clayton,  2011a;  Teodosiu  &  Sidoroff,  1976),  implying  existence 
of  a  stress  free  intermediate  or  natural  configuration  (Rajagopal  &  Srinivasa,  1998)  from  which  the  material  exhibits  an 
instantaneous  thermoelastic  response.  Typical  crystal  plasticity  models  of  high  rate  behavior  have  used  the  elastic  Green 
strain  measure  E  =  \  (FETFE  -  1)  in  their  nonlinear  elastic  stress-strain  relations  (Clayton,  2005a,  2005b;  Luscher,  Bronkhorst, 
Alleman,  &  Addessio,  2013;  Vogler  &  Clayton,  2008;  Winey  &  Gupta,  2006).  In  Section  3.1  of  the  current  work,  logarithmic 
thermoelastic  strain  e  =  lnl/E  is  used  as  a  state  variable  in  the  internal  energy,  where  FE  =  REUE  =  VERE.  Inelastic  deforma¬ 
tion  Fp  may  result  from  a  host  of  physical  mechanisms  in  single  crystals,  including  dislocation  slip  (Clayton,  McDowell,  & 
Bammann,  2004;  Teodosiu  &  Sidoroff,  1976),  deformation  twinning  (Clayton,  2009),  pore  collapse  (Barton,  Winter,  &  Reaugh, 
2009;  Clayton,  2008),  and/or  cleavage  fracture  on  preferred  planes  (Aslan,  Cordero,  Gaubert,  &  Forest,  2011;  Clayton,  2006). 
Thermodynamically  consistent  elastic-plastic  models  incorporating  logarithmic  strain  measures  have  been  developed  else¬ 
where  for  isotropic  solids  (Xiao,  Bruhns,  &  Meyers,  2007);  several  known  previous  logarithmic  models  for  anisotropic  elas¬ 
tic-plastic  crystals  (Barton  et  al„  2009;  Clayton  &  Becker,  2012)  have  posited  constitutive  equations  for  deviatoric  stress  and 
pressure  directly  (the  latter  with  an  EOS),  in  a  way  not  necessarily  consistent  with  existence  of  a  hyperelastic  total  energy 
potential. 

Solution  of  the  planar  elastic-plastic  shock  problem  is  derived  in  Section  3.2,  wherein  a  rate  independent,  but  history 
dependent,  deformation  system-level  shear  strength  model  is  applied  in  the  context  of  the  logarithmic  theory.  The  Ran- 
kine-Hugoniot  conditions  (Germain  &  Lee,  1973)  and  constitutive  relations  yield  a  set  of  coupled  nonlinear  algebraic  equa¬ 
tions  that  can  be  solved  iteratively  for  the  themomechanical  state  downstream  from  a  plastic  shock,  with  the  upstream  state 
corresponding  to  the  elastic  precursor.  The  analysis  extends  a  prior  treatment  of  isotropic  solids  (Perrin  &  Delannoy-Coutris, 
1983).  Previous  analytical  solutions  for  elastic-plastic  wave  propagation  in  crystals  have  been  restricted  to  small  strain  the¬ 
ory  (linear  elasticity)  and  isentropic  conditions  (Johnson,  1972, 1974;  Johnson,  Jones,  &  Michaels,  1970);  the  present  work,  in 
contrast,  incorporates  large  deformation,  nonlinear  themoelasticity,  and  entropy  production.  In  Section  3.3,  model  predic¬ 
tions  are  compared  with  experimental  data  on  sapphire,  diamond,  and  quartz  shocked  in  pure  mode  directions  above  the 
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HEL,  demonstrating  accuracy  of  the  model  when  proper  inelastic  deformation  mechanisms  are  incorporated  for  each  crystal 
type.  Conclusions  are  given  in  Section  4,  and  the  logarithmic  pressure-volume  EOS  is  discussed  in  the  Appendix  A. 

Notation  of  continuum  physics  is  used:  vectors  and  tensors  are  bold,  scalars  and  scalar  components  are  italic,  and  sum¬ 
mation  holds  for  repeated  subscripted  indices  (Cartesian  coordinates).  Superscripts  T,  — 1,  and  -T  denote  the  transpose,  in¬ 
verse,  and  inverse-transpose.  Scalar  and  tensor  products  obey  a  b  =  aKbK,  A  .  B  =  AyBij,  and  (a  ®  b),j  =  a/fa/. 

2.  Anisotropic  logarithmic  thermoelasticity 

2.3.  Continuum  theory 

Spatial  and  referential  coordinates  are  related  by  the  motion  x  =  x(X.  t ).  The  deformation  gradient  F  and  its  determinant 


are 

F  =  V 0x  (Fy  =  djXi),  J  =  detF  =  ^  >  0.  (2.1) 

Initial  and  current  mass  densities  are  p0  and  p.  Let  P  and  cr  denote  first  Piola-Kirchhoff  and  Cauchy  stress: 

P  =J<tF-t ■  Pj  =](J,kFJk  =^aikdkXj.  (2.2) 

Let  v  =  x  be  particle  velocity,  where  the  superposed  dot  is  a  material  time  derivative.  Balances  of  linear  and  angular  momen¬ 
tum  in  the  absence  of  body  force  are 

Vo  •  P  =  p0i>;  djPij  =  p0Xi ■  (2.3) 

PFt  =  FPt;  PijFy  =  PkjFjj.  (2.4) 

Let  ¥  and  U  denote  Helmoltz  free  energy  and  internal  energy  per  unit  initial  volume,  and  let  0  and  //  denote  absolute  tem¬ 
perature  and  entropy  density;  then 

U  =  ¥  +  6t],  (2.5) 

For  homogeneous  thermoelastic  solids,  in  general, 

¥=¥(F,6),  U  =  U  (F,  r\).  (2.6) 

Dependence  on  F  will  be  replaced  later  by  dependence  on  a  logarithmic  strain  measure  that  respects  rotational  invariance  of 
the  thermodynamic  potentials. 

The  local  balance  of  energy,  in  the  absence  of  scalar  heat  sources,  is 

U  P  F  V„  •  Q:  U  =  PijdjXt  -  djdp  (2.7) 

with  Q  heat  flux  in  reference  coordinates.  Local  entropy  production  obeys 

f\  +  V0  •  (Q/0)  ^0;  0/)  +  d/Q,  -  (Q_jdjd)/Q  >  0.  (2.8) 

Using  (2.5)  and  (2.7)  in  (2.8), 

P  :  F  -  t]9  -  ¥  -  (Q  ■  Vo0)/0  ^  0.  (2.9) 

Then  from  the  first  of  (2.6), 

(P  -  dW/0F)  :  F  -  (r)  +  d¥/d 0)0  -  Q  ■  Vo0  ^  0.  (2.10) 

from  which  the  standard  thermoelastic  constitutive  equations  of  hyperelasticity  follow: 

P  =  d¥/dF ,  i]  =  -d¥/dd.  (2.11) 


From  (2.5)  and  (2.6),  with  0  =  0(F,  /;), 

dU  d¥  d¥  <90  86  8U  8¥  86  86 

8F~8F+l)ddF  +  tldF'  dri  ~  50  8ri  +  r,dt]  +  9'  ^  ^ 

Then  from  (2.11),  in  terms  of  internal  energy,  stress  and  temperature  obey 

P  =  8U/8F,  6  =  8X8/87],  (2.13) 

Applying  the  polar  decomposition  to  the  deformation  gradient, 

F  =  RU  =  VR,  RR1  =  1,  U  =  UT,  V  =  VJ . 


(2.14) 
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Let  C  =  FJF  =  U2.  Material  logarithmic  strain  e(X,  t)  is  defined  as 

e  =  lnl/  =  llnC;  e9  =  (In  U)v  =  efl.  (2.15) 

General  definitions,  identities,  and  means  of  calculation  of  the  logarithm  of  a  second-order  tensor  are  given  in  Jog  (2008);  in 
particular, 

J  =  v/detC  =  exp  (tre)  =  exp  (e^).  (2.16) 


Free  and  internal  energies  (2.6)  are  now  posited  of  the  form 

¥=\ ¥{e,8),  U=U(e,ti).  (2.17) 

These  seem  to  have  been  rarely  used  for  anisotropic  solids,  an  exception  being  analysis  of  higher-order  moduli  in  Dluzewski 
(2000).  Conjugate  thermodynamic  variables  are 

S  =  d¥/de  =  dU/de,  rj  =  -d¥/d6,  8  =  dU/dr].  (2.18) 

From  (2.2),  the  first  of  (2.13),  (2.15),  and  the  chain  rule,  Cauchy  stress  is 


,  dU 


OF 


dU  dC 


dC  '  OF 


FT  =  2j'F§F'=j 


i  I  dU  din C 

de  '  dC 


rT  =j~1f(s  :  m)ft, 


(2.19) 


where  ¥  can  be  used  alternatively  in  place  of  LI;  in  index  notation, 

Oiy  =J  1FiKFjLMijia.Sij. 


(2.20) 


Fourth-order  tensor  M  obeys  (Jog,  2008) 


din  C  1  „  T  ^  In  A,-  -  In  A,-  _T 
1=1  1  1=1  1  J 


dC 


Here  A,-  =  2?  are  the  principal  values  of  C,  (A  B)IJKL  =  AiKBjL,  and 


P,-=  n  (C  -  Ajl)/(Aj  -  Aj). 
i=ia"w 


(2.21) 


(2.22) 


Noting  that  principal  stretches  2,-  are  eigenvalues  of  U  (and  V), 

3  1  3 

e  =  ^P,ln2i  =  ^PilnA.  (2.23) 

i=l  ^  i=l 

Let  c(e,0)  denote  specific  heat  per  unit  reference  volume  at  constant  deformation,  where  from  (2.5)  and  (2.11): 

c  =  dU/dB  =  9(dr]/d8)  =  -8(d2¥/d82).  (2.24) 

The  rate  of  internal  energy  can  be  expanded  as 

U  =  (, dU/dF )  :F  +  (dU/dri)ri  =  P  :  F+  0[(dr\/dF)  :  F  +  (dtj/d8)8\.  (2.25) 

Substituting  (2.24)  and  (2.25)  into  (2.7)  leads  to  the  temperature  rate  equation: 

c8  =  8(d2¥/dFd8)  :F-Vo  Q  (2.26) 

Defining  thermal  stress  coefficients  p(e,  8)  as 

p  =  drj/de  =  -dS/d8  =  -d2¥/ded8,  (2.27) 

(2.26)  can  be  written  as 

c8  =  -8p  :  e  —  V0  ■  Q  =  -c0f  :  c  -  V0  ■  Q.  (2.28) 

The  second-order  tensor  of  Griineisen  parameters  is 

f  =  p/c.  (2.29) 

The  following  Maxwell-type  equalities  can  be  derived  using  standard  thermodynamic  procedures  like  those  in  Thurston 
(1974)  and  Clayton  (2011a): 

8t  =  (8/c){dt]/de)  =  —dS/dri  =  -88 /de,  (2.30) 

{d/(?)7.  =  {8/c?){deld8)  =  de/dr]  =  -88/dS,  (2.31) 

«  =  de/dB  =  drj/dS.  (2.32) 
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Thermal  expansion  coefficients  are  ay.  Specific  heat  at  constant  stress  cs  obeys 

0s  =  9(drj/d6)\s  =  c+6dn:~p. 

Isothermal  and  isentropic  second-order  thermodynamic  elastic  coefficients  are 

a2  LI 
deydeia 


_  DSkl 
]>KL  ~  de„ 


f-r\  _  9Skl 

^ KL  “  det 


deydeia  aey 

Thermal  expansion  and  thermal  stress  coefficients  are  related  by 
p  =  (dr!/dS)\e :  (dS/de)  |,  =  a  :  Ce;  p,j  = 

Isentropic  and  isothermal  coefficients  are  related,  using  Maxwell  relations,  by 

c''  =  (ds/de)\0  +  (os/ae  y  ®  (ae/ae) \0  =  c"  +  ( e/c)p®p , 

or  in  indicial  notation, 

C  V  =  C;V  +  (0/c)A;fe- 


(2.33) 

(2.34) 

(2.35) 

(2.36) 

(2.37) 


Application  is  focused  on  wave  propagation  problems,  for  which  the  internal  energy  potential  proves  more  convenient 
than  the  free  energy  potential.  An  unstrained  reference  state  is  defined  by  {e,0,r\)  =  (O,0o,i]o);  temperature  and  entropy 
changes  from  this  reference  state  are  A 0  =  6  -  0o  and  A;;  =  i]  -  i]0.  Greek  subscripts  denote  Voigt  notation  for  symmetric 
indices,  e.g.,  (■)„  =  (■),,  <-►  (•)«  : 

11  <->1,  22  <->  2,  33  <->  3,  23  =  32  <-> 4,  13=31<-5,  12  =  21^6.  (2.38) 


Following  standard  convention  (Clayton,  2011a;  Thurston,  1974),  shear  strain  components  contain  a  factor  of  two,  but 
stress  and  stiffness  coefficients  do  not.  For  example,  e6  =  2ei2,S6  =  Si2,  and  C66  =  C]2i2.  Internal  energy  per  unit  reference 
volume  is  expressed  as  a  series  expansion  about  energy  U0  from  the  reference  state: 


U  —  U  o  +  Qes,  +  —  Ca£cac0 


4! 


raeaAt7  +  —TapexepA rj  -  h(rj) 


(2.39) 


Letting  (  )|0  =  (•)|e=oi?=/;o,  material  coefficients  in  (2.39)  are  constants  evaluated  at  the  reference  state,  which  is  prescribed 
stress  free: 


Uo  =  0(0,  f]0),  C.  =  (0&/aea)|o  =  O; 


C 


a/! 


'  a2u  \ 

deadei!  J 


i  - 

0 


(  fll°  ) 

Xde^depdey) 


o 


(2.40) 

(2.41) 


— 


9  dr] 
c  <9ea 


a2u 

drjdec 


BoCxp  —  I  9 


i  ara 
de,. 


cPU 


ar]dexdep 


With  constant  specific  heat  c0  =  (dU/dd)\0,  to  second  order  in  entropy  change, 

1  , 

h  =  c0[exp(A?7/c0)  -  1]  w  At/ +  (Aiy)  /c0. 


(2.42) 


(2.43) 


Thermoelastic  properties  are  usually  reported  in  the  context  of  conventional  F-based  Lagrangian  elasticity  theory 
(Clayton,  2011a;  Thurston,  1974;  Wallace,  1972),  where 

E  =  ^(C-1)  =  *(FTF-1),  U  =  U(E, ;/).  (2.44) 

Let  quantities  with  an  overbar  [i.e„  (■)]  refer  to  those  measured  with  respect  to  F-based  theory.  It  can  be  shown  (Clayton, 
2013;  Dluzewski,  2000;  Perrin  &  Delannoy-Coutris,  1983)  that  second-order  isentropic  elastic  constants  Cxp  and  Griineisen 
constants  Fa  should  be  equal  when  the  reference  state  is  unstressed  for  F-based  theory  and  e-based  (i.e.,  logarithmic) 
theory: 

C-ap  =  Ca„  ra  =  ra.  (2.45) 

This  result  is  consistent  with  the  requirement  that  UkO  when  strains  are  small.  Third-order  isentropic  constants  are 
related,  in  full  tensor  notation,  by  Dluzewski  (2000) 
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QjKLMN  —  C IJKLMN  +  2  ( J IJKLPQ  CpQMV  +  J  KLMNPQ^PQIJ  +  ^MNIJPQ^PQKl)  j 

JlJKLMN  =  g  (^IK^JM^LN  +  ^IK^JN^LM  +  ^IL^JM^KN  +  ^IL^JN^KM  +  ^lM^JK^LN  +  ^IM^JL^KN  +  ^IN^JK^LM  +  ^IN^JL^Km)- 


(2.46) 

(2.47) 


C/jkl  =  j (dCy/dEia  +  <9rja/<9£jj)|0  «  ^ (TyiS/a.  +  C/aSy), 


A  standard  assumption  (Wallace,  1980)  for  weak  shocks  is  pr  rj  constant,  which  yields  (Clayton,  2013) 

(2.48) 

-'ijkl  are  then  related  using  the  procedure  outlined  in  Clayton  (2013),  giving 

(2.49) 


Higher-order  Griineisen  parameters  f 


a/J  ' 


Tyia  =  r,JKL  +  ^(r,Kd]L  +  TjKdiL  +  r,LSjK  +  r]Ls,K) »  ^(TySKL  +  rKLSy  +  rIKsJL  +  r]KsIL  +  r ILdjK  +  rJLSiK). 


2.2.  Planar  shock  solution 

A  shock  wave  is  represented  as  a  propagating  surface  across  which  there  may  exist  jump  discontinuities  in  mass  density, 
particle  velocity,  strain,  stress,  entropy,  temperature,  and  internal  energy.  Considered  here  are  1-D  (i.e,  normal,  planar,  or 
longitudinal)  shocks.  Quantities  associated  with  material  ahead  of  the  shock  are  denoted  with  superscript  +,  with  material 
behind  superscript  In  the  present  analysis  of  Section  2.2,  material  ahead  of  the  shock  is  at  rest,  undeformed,  unstressed, 
and  at  ambient  reference  temperature  00;  these  assumptions  are  removed  later  in  Section  3.2.  The  jump  in  an  arbitrary  quan¬ 
tity  (■)  across  the  shock  is  written 

[(•)]  =  (-r-(0+-  (2.50) 

The  shock  moves  at  steady  natural  velocity  D  >  0  in  the  X  =  Xi  direction.  The  deformation  gradient  is 

"F  0  0]  ri+{  0  0" 

F  =  [F,y]“  =010  =  0  1  0  ;  F+  =  1.  (2.51) 

_0  0  lj  [  0  01 

Behind  the  shock,  with  x  =  x^  and  u  =  if  particle  coordinate  and  displacement, 

F=§=1+S  =  1  +  ^-r =  v~0  =  r’  i  =  du/dX ■  (2-52) 

Attention  is  restricted  to  compressive  shocks,  for  which  0  <  F  sg  1  and  —1  <  £  ^  0.  The  only  nonzero  component  of  logarith¬ 
mic  strain  e  is 


e  =  eu  =  InF  =  ln(l  +  f). 


(2.53) 


The  “shock  stress”  or  “shock  pressure”  is  the  longitudinal  force  per  unit  reference  area  (or  equivalently,  current  area)  behind 
the  shock,  positive  in  compression: 


P  =  -Pu  =  -KFn°rk)  =-au- 


(2.54) 


Let  p  =  p~  and  0=07  be  mass  density  and  particle  velocity  in  the  shocked  state.  Conservation  laws  for  mass,  linear 
momentum,  and  energy  -  i.e.,  particular  forms  of  the  Rankine-Hugoniot  equations  -  are,  respectively  (Clayton,  2013; 
McQueen  et  al„  1970;  Thurston,  1974), 


p0t>  =  p(D  ~v)<=^i=  -d/d, 
P  =  p0T>V 

n 


Pu  =  VI-p0v 


Po&  =  -P/( 

2 1  m 


p0v2  =  -P{, 

1 


(2.55) 

(2.56) 

(2.57) 


From  (2.55),  requiring  1  >J>  0  leads  to  constraints  D  >  d  >  0.  In  (2.57),  the  usual  adiabatic  assumption  of  null  heat  con¬ 
duction  has  been  used.  Shock  compression  is  neither  isothermal  nor  isentropic;  the  entropy  inequality  is  (Germain  &  Lee, 
1973) 


b /Pol  >  0 


>  0. 


(2.58) 


Subsequent  derivations  invoke  internal  energy  U  =  U(e,  if)  of  (2.39).  Derivatives  of  U  with  respect  to  strain  depend  only 
on  entropy  changes  A; 7  from  the  reference  state  and  are  independent  of  ;/0  =  r/+,  and  stress  and  temperature  depend  only  on 
derivatives  of  internal  energy  with  respect  to  strain  and  entropy  and  are  independent  of  U0.  Therefore,  let 

U0  =  U+  =  0,  ;/0  =  rj+  =  0  =►  [U]  =  U  =  U,  M  =  tf  =  Ai /  =  rj-  (2.59) 

0+  =  {dU/dr])+  =  0O,  8-  =  (, dU/dri )~  =  8  =>  [0]  =  A0.  (2.60) 
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Stress  components  thermodynamically  conjugate  to  e  are  related  to  P  via  (2.19),  (2.21 ),  (2.22)  and  (2.54),  which  reduce  here 
to 


P  =  ~J  'FvFwSklMkljn  =  -FSiiMim  =  -S/(l  +  £), 

(2.61) 

S  =  S„  =dU/deu  =dU/de, 

(2.62) 

where  all  quantites  are  evaluated  behind  the  shock. 

Expanding  {  and  (1  +  f)_1  in  series  to  order  five  in  logarithmic  strain  e  =  eu, 

i  =  exp  (e)  -  1  «e  +  ^e2  +  ^e3+ie4+ie5, 

(2.63) 

1/(1 +  ^“1  -e  +  le2-le3+le4-le5, 

(2.64) 

^/(l+0-e-^2  +  ^e3-^e4  +  T^e5. 

(2.65) 

Using  (2.56),  (2.59),  (2.61),  and  (2.65),  balance  (2.57)  becomes,  in  terms  of  e, 

u = -  ^ + «  -  lie  -  + le3  4e4 + iiie5)  • 

(2.66) 

Internal  energy  function  (2.39)  -  using  (2.43)  and  specialized  to  uniaxial  strain  with  (2.53)  and  (2.59),  to  first  order  in  en¬ 
tropy  -  becomes 


1 


1 


1 


=  2  Cue  +  g Cme  +  —  Cmie  -  0O(  r,ei7  +  -rne^7  -  t] 


1 


(2.67) 


When  0  is  a  linear  function  of  entropy  as  in  (2.67),  then  a  solution  for  i](e)  can  be  obtained  analytically  in  closed  form,  as  in 
what  follows.  This  form  of  internal  energy  is  most  valid  for  weak  to  moderate  elastic  shocks  wherein  i)/c0  <c  2,  which  will  be 
verified  a  posteriori  in  later  examples.  Conjugate  thermodynamic  stress  is 


S  —  dU / 9e  —  Cn  &  +  ^  CmC2  +  c  61111C3  —  dojTi  +  rne)fj- 


(2.68) 


2  6 

This  is  substituted  into  (2.66),  which  then,  when  considered  with  (2.67),  provides  two  equations  in  three  unknowns  ( U ,  e,  ij). 
Writing  17(e)  as  a  polynomial  with  constant  coefficients  a0,  a  1 ,  a2, . . ., 


i]  =  a0  +  cq  e  +  a2e2  +  a3e3  +  a4e4  +  a5e5  +  •  ■  • 


(2.69) 

Substituting  (2.69)  into  (2.66)  and  (2.67),  setting  U  =  U,  equating  coefficients  of  like  powers  of  e  up  to  order  five,  and  noting 
fj0  =  i)( 0)  =  0  from  (2.59), 


do  =  =  <h  =  0,  <33  =  ®o'(— 3Cn  +  Cm), 

04  =  2^S0^[12C^]  —  6Cm  +  Cmi  +  r,(-3C„  +  Cm)], 

Q5  =  Hq  1  [  Ci  1  +  2 Ci  1 1  —  2Ci in  +  r,(9Ci,  —  5Ci  1 1  +  Cmi)  +  r2(— 3Cn  +  Cm)]. 
Substitution  of  entropy  r\(e),  now  known  to  fifth  order  in  strain,  into  (2.68)  gives 

S(e)  =  Cue  +  -  Cme2  -(-  ^-Cim  —  doria3J e3  —  0o(riQ4  +  rna3)e4  —  0o(riQ5  +  rn04)e5. 


(2.70) 

(2.71) 

(2.72) 

(2.73) 


Use  of  this  result  for  stress  with  (2.61)  and  Hugoniot  equations  (2.55)-(2.57)  then  gives  shock  pressure,  internal  energy,  par¬ 
ticle  velocity,  and  shock  velocity  in  terms  of  single  variable  e  =  In)  =  ln^: 


w 

P=  -S/exp  (e).  [U]  =  ^S[l  -  exp  (— e)], 

^  =  {(-S/Po)!1  -  exP  C-e)]}1/2, 

£  =  {(VPo)[l  -  exp  (-e)]}’/2[l  -  exp  (e)]. 

Finally,  thermoelastic  temperature  rise  due  to  entropy  production  is  simply 

9  =  dU/dr)  =  -  r^-^Tiie2 


(2.74) 

(2.75) 

(2.76) 

(2.77) 


In  the  limit  of  very  low  shock  stress,  shock  velocity  approaches  the  longitudinal  linear  elastic  wave  speed 

Co  =  (C„/p0),/2  (2.78) 
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and  shock  stress  in  the  first  of  (2.74)  approaches  the  isentrope 


P'1  =  -exp(-e)^Cne  +  iciiie2  +  gCime3^. 


(2.79) 


2.3.  Results:  elastic  solution 

Theory  and  analytical  solutions  derived  in  Section  2.2  are  now  applied  to  analyze  shock  compression  behavior  of  single 
crystals  of  three  hard  minerals:  sapphire  (a-Al203  or  corundum),  diamond  (C),  and  quartz  (oc-Si02).  These  materials  are 
considered  because  elastic  deformations  in  excess  of  several  percent  volumetric  compression  can  be  achieved  in  uniaxial 
compression  prior  to  any  inelastic  deformation  that  could  render  the  analysis  of  Section  2.2  unrealistic.  These  materials  also 
belong  to  the  limited  set  of  anisotropic  crystals  whose  complete  third-order,  and  in  some  cases  fourth-order,  elastic  con¬ 
stants  have  been  reported. 

Analytical  solutions  are  studied  for  uniaxial  shock  compression  involving  internal  energy  function  (2.39)  incorporating 
logarithmic  strain  e  with  elastic  constants  up  to  possibly  fourth  order.  Certain  results  are  also  compared  with  those  obtained 
previously  (Clayton,  2013)  for  Lagrangian  (E-based)  and  Eulerian  (D-based)  nonlinear  thermoelastic  models.  Sapphire  and 
quartz  have  trigonal  (i.e.,  rhombohedral)  symmetry  and  are  analyzed  for  compression  along  two  pure  mode  directions: 
the  a-axis  (X-cut,  [1210])  and  c-axis  (Z-cut,  [0001]).  The  b-axis  is  not  a  pure  mode  direction  in  trigonal  crystals;  the  plane 
wave  approximation  used  for  analysis  of  Y-cut  quartz  in  Clayton  (2013),  which  omits  possible  transverse  displacements, 
is  not  repeated  herein.  Diamond  is  cubic  and  is  analyzed  for  pure  mode  compression  along  a  cube  axis  (X-cut,  [100]).  Elastic 
constants  are  interchanged  as  needed  for  consistency  with  notation  of  Section  2.2.  For  example,  for  c-axis  (i.e.,  Z-cut)  com¬ 
pression,  Cn  is  replaced  by  C33,  Cm  by  C333,  Ti  by  r3,  etc. 

Thermoelastic  properties  are  listed  in  Table  1.  Second-order  elastic  constants  are  reported  from  ultrasonic  experiments 
(McSkimin,  Andreatch,  &  Glynn,  1972;  McSkimin,  Andreatch,  &  Thurston,  1965;  Winey,  Gupta,  &  Hare,  2001).  Higher-order 
constants  are  converted  from  reported  Lagrangian  values  (Hankey  &  Schuele,  1970;  Nielsen,  1986;  Thurston  et  al„  1966)  to 
logarithmic  values  via  (2.46)  and  (2.49).  Complete  sets  of  anisotropic  second-  and  third-order  constants  are  listed  for  later 


Table  1 

Thermoelastic  single  crystal  properties  (0O  =  295  I<;  Cy,s  in  GPa). 


Property 

Sapphire 

Diamond 

Quartz 

Cn 

498 

1079 

88 

Cl2 

163 

124 

7 

Cl3 

117 

C12 

12 

C14 

23 

- 

-18 

C33 

502 

Cn 

106 

C44 

147 

578 

58 

Cm 

-792 

174 

315 

C112 

-764 

-552 

-331 

C113 

-729 

C112 

36 

C114 

-9 

- 

-199 

Cl  23 

-289 

0 

-294 

Cl  24 

62 

- 

-33 

Cl33 

-688 

C112 

-288 

Cl34 

154 

- 

-16 

Cl  44 

-162 

124 

-125 

Cl  55 

-559 

-843 

-34 

C222 

-1532 

C111 

193 

C333 

-328 

Cm 

-181 

C344 

-487 

Cl55 

65 

C444 

-16 

- 

-249 

C456 

2  (Cl55  —  C144) 

-433 

2  (C155  -  C144) 

C1111 

- 

- 

to4 

C3333 

- 

- 

to4 

r, 

1.29 

0.81 

0.74 

r3 

1.29 

r. 

0.58 

r„ 

3.87 

2.43 

2.22 

r,2 

1.29 

0.81 

0.74 

r,3 

1.29 

r,2 

0.66 

r33 

3.87 

r,, 

1.74 

Bo 

254 

442 

38 

Go 

166 

538 

48 

Bo 

4.3 

4.0 

6.3 

Po  [g/cm3] 

3.98 

3.51 

2.65 

Co  [MPa/K] 

3.10 

1.73 

1.95 

Ph/Ci, 

0.05 

0.08 

0.10 
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use  in  Section  3.3;  in  fact,  only  longitudinal  constants  are  needed  presently  for  analysis  of  the  axial  purely  thermoelastic  re¬ 
sponse,  as  is  clear  from  the  solution  derived  in  Section  2.2.  Bulk  modulus  B0  and  its  pressure  derivative  BJ,  (Gieske  &  Barsch, 
1968;  McSkimin  &  Andreatch,  1972;  McSkimin  et  al„  1965)  are  also  listed;  these  enter  the  hydrostatic  EOS  described  in  the 
Appendix  A.  For  crystals  with  cubic  symmetry, 

B0  =  i(C„+2C12),  (2.80) 


while  for  crystals  with  trigonal  or  hexagonal  symmetry, 

Bo  =  [(Ci,  +  C12)C33  -  2Cy /[C„  +  C,2  +  2C33  -  4C13],  (2.81) 

Constant  B'0  is  formally  related  to  combinations  of  second-  and  third-order  elastic  constants  (Guinan  &  Steinberg,  1974),  but 
is  often  measured  directly  or  fit  to  high  pressure  data.  Also  listed  for  each  crystal  is  the  Voigt-average  shear  modulus  G0,  used 
later  in  Section  3  for  normalization  of  shear  strength.  For  crystals  with  cubic  symmetry, 

Go  =  g  (Cn  —  C12  +  3C44),  (2.82) 


while  for  crystals  with  trigonal  or  hexagonal  symmetry, 


Go  =  |  [2Cn  +  C33]  -  [Ci2  -I-  2C,3]  +  3 


2  (C11  —  Ci2)  +  2C44 


(2.83) 


As  discussed  later,  fourth-order  constants  Cmi  and  C3333  for  quartz  are  fit  to  shock  velocity  versus  particle  velocity  data 
(Fowles,  1967),  but  these  are  not  needed  for  sapphire  or  diamond  described  with  logarithmic  elastic  theory.  Maximum 
HEL  stresses  PH  from  shock  experiments  (Graham,  1972;  Fowles,  1967;  Lang  &  Gupta,  2010;  Wackerle,  1962)  are  shown 
for  reference,  normalized  by  second-order  moduli  (X-cut  shown).  The  domain  of  validity  of  elastic  analysis  can  be  estimated 
as  TL  >  Cll(rFHir  In  sapphire,  the  elastic  domain  for  Z-cut  quartz  is  about  the  same  as  that  for  X-cut  (Graham,  1972;  Graham  & 
Brooks,  1971);  in  quartz,  the  elastic  domain  for  Z-cut  is  demarcated  by  ^33  «  0.15  (Fowles,  1967;  Wackerle,  1962). 

Predicted  shock  velocity  35  versus  particle  velocity  0  is  compared  with  experimental  shock  compression  data  of  Graham 
and  Brooks  (1971 )  in  Fig.  la  for  X-  and  Z-cut  sapphire.  Experimental  data  are  obtained  from  flyer-plate  and  explosive  loading 
configurations;  in  the  latter,  two-wave  structures  arose  (Graham  &  Brooks,  1971).  Data  considered  here  correspond  only  to 
the  elastic  shock,  with  the  secondary,  slower  "plastic”  wave  in  which  the  HEL  was  exceeded  not  addressed  here;  the  plastic 
response  will  be  modeled  in  Section  3.  Velocities  are  normalized  by  wave  speed  (2.78).  Predictions  marked  “3rd  order”  are 
obtained  using  complete  solutions  and  all  material  constants  in  Table  1.  Predictions  marked  “2nd  order”  assume  Cm  =  0. 
Note  that  fourth-order  constant  Cun  is  not  needed  for  either  orientation  (Table  1 ).  Considering  scatter  in  the  data,  3rd  order 
and  2nd  order  fits  are  adequate  for  each  orientation,  giving  nearly  linear  35-t>  curves  over  the  compression  range  for  which 
sapphire  remains  elastic  (^  >  0.95).  Second  order  predictions  for  X-  and  Z-cut  orientations  are  almost  identical  and  corre¬ 
spond  to  the  broken  line  with  lowest  slope  in  Fig.  1(a).  Hugoniot  stress  (i.e„  P)  normalized  by  Cn  or  C33  is  shown  in  Fig.  1(b) 
and  (c),  along  with  experimental  data  (Graham  &  Brooks,  1971 ).  Also  shown  are  2nd  order  predictions  made  using  analytical 
solutions  given  in  Clayton  (2013)  for  Lagrangian  E-based  theory  and  Eulerian  e-based  theory.  For  each  orientation,  2nd  and 
3rd  order  logarithmic  model  predictions  provide  close  agreement  with  experiment.  The  other  nonlinear  elastic  models 
(Clayton,  2013)  are  comparatively  inaccurate,  with  2nd  order  Eulerian  theory  too  stiff  and  2nd  order  Lagrangian  theory 
too  compliant. 

Predictions  of  normalized  shock  velocity  and  Hugoniot  stress  in  diamond  are  given  in  Fig.  2(a)  and  (b),  respectively,  com¬ 
pared  with  experimental  data  of  Lang  and  Gupta  (2010).  Fourth-order  constant  Cun  is  not  needed.  From  Fig.  2(a),  differences 
in  2nd  and  3rd  order  predictions  of  logarithmic  theory  are  small,  with  both  providing  close  agreement  with  experimental 
velocity  data.  Likewise  in  Fig.  2(b),  differences  in  2nd  and  3rd  order  predictions  of  logarithmic  theory  are  small,  with  both 
providing  close  agreement  with  experimental  stress  data.  The  other  2nd  order  elastic  models  (Clayton,  2013)  do  not  accu¬ 
rately  predict  Hugoniot  stress,  with  2nd  order  Eulerian  theory  too  stiff  and  2nd  order  Lagrangian  theory  too  compliant. 

Predicted  shock  velocity  35  versus  particle  velocity  u  is  compared  with  experimental  shock  compression  data  of  Fowles 
(1967)  in  Fig.  3a  for  X-  and  Z-cut  quartz.  Experimental  data  are  obtained  from  plane-wave  explosive  loading  (Fowles, 
1967);  data  considered  here  correspond  only  to  the  first,  elastic  shock  in  each  test.  Fourth-order  constant  Cun  was  fit  to 
the  data  in  Fowles  (1967).  Predictions  marked  “3rd  order”  assume  Cmi  =  0.  Predictions  of  shock  stress  are  compared  with 
those  of  3rd  order  Lagrangian  and  Eulerian  nonlinear  elastic  solutions  (Clayton,  2013)  in  Fig.  3(b).  For  each  orientation,  4th 
order  theories  are  required  to  most  accurately  match  the  experimental  Hugoniot  data;  2nd  and  3rd  order  models  are  all  too 
compliant.  However,  3rd  order  logarithmic  theory  more  closely  matches  the  data  for  Z-cut  quartz  than  3rd  order  Lagrangian 
and  Eulerian  theories. 

Predictions  of  the  present  logarithmic  nonlinear  thermoelasticity  theory  for  temperature  rise  0  (normalized  by  reference 
temperature  60),  entropy  jump  across  the  shock  r/  (normalized  by  specific  heat  c0),  and  Hugoniot  stress  P(normalized  by  uni¬ 
axial  isentropic  stress  P’1)  are  listed  in  Table  2.  Temperatures  are  computed  from  (2.77),  entropy  from  (2.69),  and  isentropes 
from  (2.79),  for  each  crystal  with  the  full  set  of  material  properties  given  in  Table  1.  Predicted  temperature  rise  is  fairly  small 
for  elastic  shock  loading,  similar  to  results  reported  for  Lagrangian  and  Eulerian  theories  in  Clayton  (2013).  Entropy  is  posi¬ 
tive  in  agreement  with  (2.58)  and  is  of  the  same  order  of  magnitude  reported  for  Lagrangian  and  Eulerian  theories  in  Clayton 
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u/C0 

(a) 


VIV0  VIV0 

(b)  (c) 


Fig.  1.  Analytical  nonlinear  elastic  solutions  and  experimental  data  (Graham  &  Brooks,  1971 )  for  sapphire:  (a)  shock  velocity  vs.  particle  velocity,  X-  and  Z- 
cut  (b)  shock  stress  vs.  volume  ratio,  X-cut  (c)  shock  stress  vs.  volume  ratio,  Z-cut. 


o° 


u/C0  VIV0 

(a)  (b) 


Fig.  2.  Analytical  nonlinear  elastic  solutions  and  experimental  data  (Lang  &  Gupta,  2010)  for  diamond:  (a)  shock  velocity  vs.  particle  velocity,  X-cut  (b) 
shock  stress  vs.  volume  ratio,  X-cut. 
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(a) 


(b)  (c) 


Fig.  3.  Analytical  nonlinear  elastic  solutions  and  experimental  data  (Fowles,  1967)  for  quartz:  (a)  shock  velocity  vs.  particle  velocity,  X-  and  Z-cut  (b)  shock 
stress  vs.  volume  ratio,  X-cut  (c)  shock  stress  vs.  volume  ratio,  Z-cut. 


Table  2 

Thermodynamic  predictions  of  nonlinear  logarithmic  theory. 


Material 

Direction 

V/V0 

e/e o 

n/c0 

p/p-i 

Sapphire 

X 

0.96 

1.049 

0.016 

1.0008 

0.92 

1.094 

0.145 

1.0028 

0.88 

1.133 

0.556 

1.0056 

Z 

0.96 

1.049 

0.012 

1.0006 

0.92 

1.094 

0.110 

1.0022 

0.88 

1.133 

0.420 

1.0045 

Diamond 

X 

0.96 

1.031 

0.036 

1.0003 

0.92 

1.059 

0.327 

1.0011 

0.88 

1.083 

1.248 

1.0023 

Quartz 

X 

0.96 

1.028 

0.001 

1.0002 

0.92 

1.054 

0.030 

1.0013 

0.88 

1.076 

0.178 

1.0041 

Z 

0.96 

1.022 

0.007 

1.0005 

0.92 

1.042 

0.086 

1.0021 

0.88 

1.060 

0.400 

1.0046 

(2013).  Recall  from  Section  2.2  that  the  present  analytical  solution  assumes  the  contribution  to  internal  energy  from  entropy 
is  linear  in  (2.67),  most  accurate  for  ;/  <  2 c0.  From  Table  2,  such  conditions  hold  for  ^  >  0.92.  Examination  of  stresses  in 
Table  2  shows  that  ^  <  1.01  in  all  cases,  meaning  an  isentropic  assumption  for  elastic  shock  compression  is  accurate  for  pre¬ 
dicting  axial  stress  up  to  the  HEL  in  these  crystals.  Values  listed  in  Table  2  are  considered  extrapolations  when  compression 
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exceeds  the  HEL.  Above  the  HEL,  a  nonlinear  thermoelastic-plastic  theory  incorporating  inelastic  deformation  mechanisms 
and  associated  dissipation,  as  developed  and  applied  next  in  Section  3,  becomes  necessary. 

3.  Anisotropic  logarithmic  thermoelastic-plastic  theory 

3.3.  Continuum  theory 

Assuming  again  that  x  is  differentiable,  deformation  gradient  (2.1)  is  decomposed  multiplicatively  (Clayton,  2011a; 
Teodosiu  &  Sidoroff,  1976)  as 

F=V0x  =  FEFp;  FiJ  =  djXi  =  FfKFpK].  (3.1) 

Here  FE  and  Fp  denote  deformation  “gradient”  mappings  associated  with  thermoelasticity  and  defects  (e.g.,  dislocation  slip, 
twinning,  or  cleavage  fracture),  though  neither  generally  anholonomic  mapping  (Clayton,  2012)  need  be  compatible  (i.e„  a 
true  gradient  of  a  vector  field).  Each  is,  however,  presumed  to  have  positive  determinant:  JE  =  detFE  >  0  and  f  =  detFp  >  0. 
The  usual  stress  definitions  and  local  balance  laws  of  continuum  mechanics  of  Section  2.1  still  apply  [(2.2),  (2.3),  (2.4),  (2.7), 
and  (2.8)],  again  here  assuming  no  body  forces  or  scalar  heat  sources. 

Logarithmic  thermoelastic  strain  is  now  defined  as  (no  superscript) 

e  =  lnl/E  =  ^lnCE,  (3.2) 

where  the  polar  decomposition  FE  =  REL3E,  and  CE=FETFE.  Analogously  to  (2.16),  thermoelastic  volume  change  obeys 
/  =  exp  (tre). 

Let  c  represent  a  generic  internal  state  variable  associated  with  evolution  of  microstructure,  for  example,  dislocation  or 
crack  density.  Here,  f  is  assumed  scalar,  but  generalization  to  higher-order  tensors  (Bammann,  1984)  poses  no  theoretical 
difficulties.  Assuming  uniform  properties  in  the  reference  state,  free  and  internal  energy  densities  per  unit  initial  volume 
are  of  the  form 

y=  V(e,e,i),  U=0(e,  »/,£),  (3.3) 

related  as  usual  via  U  =  lP  +  Otj.  Define  thermodynamic  conjugate  forces  as 

S  =  dU/de,  f  =  -d(J/dt  (3.4) 

Using  (2.7),  (2.8),  (3.3),  and  (3.4),  and  considering  admissible  thermomechanical  processes,  the  following  constitutive  laws 
can  be  derived  consistently  with  the  first  and  second  laws  of  thermodynamics  (see  e.g.,  Clayton,  2011a): 

<t=/-1Fe(S  :  M)FET,  i\  =  -dW/dd,  6=dU/dr].  (3.5) 

Here,  M  =  d\n CE/9CE  is  computed  analogously  to  (2.21).  With  Lp  the  inelastic  velocity  gradient,  entropy  inequality  (2.8)  re¬ 
duces  to 

Jo  :  (FeLpFe-’)  +  -  Q  ■  Vo0  >  0;  (Lp  =  FpFm). 

Defining  specific  heat  at  constant  elastic  strain  as  in  (2.24),  i.e.,  c  =  dU/86 

c8  =  Jo (FeLpFe-j)  +  8(dS/dB)  :  e  +  [c  -  8{dt/d0)}t  -  V0  •  Q. 

Generic  kinetic  equations  for  inelasticity  are  of  the  state-dependent  form 

Lp  =  Lp(Fe,  ij,  {),  £  =  i(FE,»7,a  (3.8) 

Thermodynamic  definitions  and  identities  in  (2.27)  and  2.29,  (2.30)-(2.36)  and  (2.37)  still  apply  in  the  elastic-plastic  case, 
where  now  it  is  understood  that  e  is  the  thermoelastic  part  of  the  total  strain  as  in  (3.2). 

For  anisotropic  single  crystal  plasticity,  contributions  to  inelastic  deformation  can  often  be  attributed  to  deformation  on 

distinct  planes  with  reference  unit  normal  vector  m$,  in  unit  reference  direction  sg,  where  oc  =  1,2, _ n,  with  n  the  number 

of  possible  deformation  systems.  Let  ;f  denote  the  scalar  deformation  rate  on  system  a,  which  in  general  here  can  be  due  to 
dislocation  glide,  deformation  twinning,  or  cleavage  fracture.  The  inelastic  velocity  gradient  and  inelastic  dilatation  rate 
become 

lp  =  E^so 0  m«  ip=/E^so-mo-  (3-9) 

a  a 

Lattice  directors  deform  thermoelastically  via 

sa  =  FeSq,  ma  =  FE~TmrQ.  (3.10) 

Let  =  / yadt  denote  the  cumulative  deformation  on  system  a  over  a  time  increment  in  which  deformation  increases 
monotonically.  For  dislocation  glide,  then  represents  isochoric  plastic  shear  (sg  ■  mg  =  0);  for  twinning,  yx  =  yT/a,  with 


(3.6) 

=  -0d2'F/8d2,  the  balance  of  energy  becomes 

(3.7) 
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/“  and  yT  the  local  volume  fraction  of  twinned  crystal  and  characteristic  twinning  shear  (Clayton,  2009,  2010c)  (twinning  is 
likewise  isochoric  with  sjj  •  mg  =  0).  For  cleavage  (Clayton,  2010b),  mode  II/I1I  fracture  is  isochoric  (sg  ■  mg  =  0),  but  mode  I  is 
not  (sg  ■  mg  =  1  and  ya  3:  0.)  Dislocation  glide  is  lattice  preserving  and  does  not  affect  thermoeiastic  properties  such  as  elas¬ 
tic  moduli.  Deformation  twinning  induces  lattice  rotation  that  can  affect  moduli  and  slip  directors  within  twinned  regions. 
Fractures  may  induce  degradation  of  thermoeiastic  properties. 

Using  (3.9)  and  (3.10),  dissipation  (3.6)  and  temperature  rate  (3.7)  become 

+  Cf  -  0.  ■  Vo0  >  0,  (Ta=Jo  :sx®ma)-,  (3.11) 

a 

cB  =  J2r T  +  S(dS/d8)  :  e+  [C  —  B(dC/dB)]t  -  V0  Q  (3.12) 


The  thermodynamic  driving  force  for  inelastic  deformation  on  system  a  is  resolved  Kirchhoff  stress  r“.  For  rate  dependent 
inelasticity,  a  more  specific  form  of  the  first  of  (3.8)  is 

f  =  Ta,e,£).  (3.13) 

In  the  rate  independent  limit,  often  thought  adequate  for  describing  strong  crystals  such  as  ceramics  whose  shearing  resis¬ 
tance  g"  is  an  appreciable  fraction  of  the  theoretical  strength  (Graham  &  Brooks,  1971),  (3.13)  is  replaced  with 

r  a<g'z(0,{)  <=>  f  =  0,  T“  =  ga(0,£)  If  I  >0.  (3.14) 


Rate  independent  theory  also  proves  useful  as  a  limiting  case  for  describing  the  response  in  problems  such  as  shock  com¬ 
pression  wherein  the  true  strain  rate  is  unspecified  or  unknown  (e.g.,  a  shock  represented  mathematically  as  a  moving  sur¬ 
face  of  discontinuity,  as  in  Germain  &  Lee  (1973),  Perrin  &  Delannoy-Coutris  (1983)  and  in  later  Section  3.2  of  the  present 
work). 

Application  is  again  focused  on  wave  propagation  problems.  Internal  energy  per  unit  reference  volume  is  expressed  as  a 
series  expansion  about  energy  U0  from  the  reference  state  as  in  (2.39)-(2.42)  and  (2.43): 


U  —  Uo  +  +  gCaft> e«e/Jey  ~^24^'01  Pyze<xelieye6  +  9oAl] 


1 T  At]  \~yOy  r^eae^ 
ZCq  / 


+  2K<^oi2  ■ 


(3.15) 


The  final  term  represents  stored  energy  of  lattice  defects  for  example,  Regueiro,  Bammann,  Marin,  and  Garikipati  (2002)  and 
Clayton  (2005a),  with  k  =  i  <0  a  dimensionless  constant.  Here,  isentropic  moduli  are  generally  anisotropic  but  constant  as 
in  Section  2.1 ;  for  simplicity,  possible  effects  of  twinning  or  damage  on  elastic  moduli  are  omitted.  Identities  in  (2.45)-(2.48) 
and  (2.49)  still  apply. 


3.2.  Planar  shock  solution 


Consider  a  continuous  and  initially  homogeneous  slab  of  material  through  which  a  planar  shock  moves,  in  the  Xi -direc¬ 
tion,  with  natural  velocity  £.  As  in  Section  2.2,  let  superscripts  +  and  -  label  quantities  in  the  material  ahead  (i.e.,  upstream) 
and  behind  (i.e.,  downstream)  from  the  shock.  Let  [(■)]  =  (•)“  -  (-)+  and  ((■))  =  ![(•)”  +  (  )+]  denote  the  jump  and  average  of  a 
quantity  across  the  shock.  Let  n  be  a  unit  normal  vector  to  the  planar  shock,  i.e.,  n  =  dx/dx\.  The  only  nonvanishing  compo¬ 
nent  of  particle  velocity  is  u  =  v  ■  n.  The  Cauchy  stress  component  normal  to  the  shock  front  is  a  =  a  :  (n  ®  n)  =  Cn.  The  rel¬ 
ative  velocity  of  the  material  with  respect  to  the  shock  is  v  =  d  -  2>.  Let  u  =  U/p0  denote  internal  energy  per  unit  mass. 
Appropriate  forms  of  the  Rankine-Hugoniot  conditions  for  conservation  of  mass,  momentum,  and  energy  are  (Germain  & 
Lee,  1973) 


o' 

II 

3. 

(3.16) 

M  -  Pv[v]  =  0, 

(3.17) 

[pvfli  +  jv2)  -  <Tv]  =  0. 

(3.18) 

The  material  need  not  be  deformed  uniaxially  according  to  these  conditions,  but  shock  velocity  and  particle  velocity  must 
both  be  rectilinear  in  the  -direction  so  that  only  normal  traction  is  discontinuous.  Therefore,  these  equations  can  apply 
for  shock(s)  passing  through  a  pre-stressed  material  such  as  a  plastic  wave  following  an  elastic  precursor,  making  them  more 
general  than  2.55,  2.56  and  2.57  of  Section  2.2  that  apply  only  for  a  single  shock  passing  through  an  initially  unstressed  slab. 
Adiabatic  conditions  have  been  assumed  (Germain  &  Lee,  1973):  Q.  =  0,  leading  to  entropy  production  requirement  (2.58). 
Using  (3.16)  and  (3.17),  energy  conservation  condition  (3.18)  can  be  rewritten  as  (Germain  and  Lee,  1973) 

M  =  (tf}[l/p]  <=*  M  =  (ff>[/]-  (3.19) 

Assume  that  the  upstream  state  and  shock  velocity  are  known.  The  downstream  state  is  defined  by  variables  (tr,  er~.  u~). 
The  Rankine-Hugoniot  conditions  provide  three  equations  for  determining  this  state;  in  order  to  fully  obtain  the  down¬ 
stream  state,  a  fourth  equation  is  supplied  by  the  constitutive  model,  or  another  downstream  variable  must  be  known. 
For  example,  in  experiments  particle  velocity  tr  is  often  measured  in  addition  to  D. 
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Considered  in  what  follows  are  longitudinal  elastic-plastic  shocks  corresponding  to  planar  impact  in  pure  mode  direc¬ 
tions  in  single  crystals  (i.e.,  directions  parallel  to  an  axis  of  threefold  or  greater  symmetry).  A  sample  of  material  subjected 
to  a  step  or  ramp  loading  in  normal  stress  along  this  direction,  with  no  applied  shear  stress,  develops  a  two-wave  structure 
consisting  of  a  single  longitudinal  elastic  wave  (i.e.,  the  elastic  precursor),  followed  by  a  single  longitudinal  plastic  wave  if 
the  HEL  is  exceeded.  Implicitly,  overdriven  shocks  in  which  the  plastic  wave  overtakes  the  elastic  precursor  are  not  consid¬ 
ered.  Total  deformation  is  thus 


’n 

+ 

o 

o 

O 

o 

1 

Uh 

II 

0  1  0 

°  0  1 

II 

1  + 

*Ti 

1 

II 

0  1  0 

0  0  1 

(3.20) 


where  hereafter  superscripts  +  and  -  label  quantities  in  front  of  and  behind  the  plastic  shock,  i.e.,  at  the  HEL  state  of  the 
elastic  precursor  and  the  final  inelastic  state.  Note  that  at  the  HEL  state,  F+  =  (FE)  ,  while  behind  the  shock,  F  =  (FEFP)  . 
For  crystal  orientations  of  lower  symmetry,  or  if  shear  stress  is  applied  in  addition  of  longitudinal  pressure,  quasi-longitu¬ 
dinal  and  multiple  quasi-transverse  elastic  and  plastic  waves  would  appear  (Johnson,  1972);  such  more  complicated  prob¬ 
lems,  for  which  (3.20)  no  longer  applies,  are  not  addressed  herein. 

For  the  highly  symmetric  orientations  considered  in  detail  for  sapphire  and  diamond  in  Section  3.3,  n  inelastic  shear  sys¬ 
tems  are  active  simultaneously  at  shock  stresses  exceeding  PH,  all  at  the  same  rate  y  =  ya.  For  monotonic  loading,  integration 
of  (3.9)  yields  the  plastic  deformation,  to  order  three  in  y: 


Fp(y)  =  exp 


)’J2S  o  1 


m, 


(3.21) 


with  cumulative  slip  y  to  be  determined  in  the  analysis.  This  slip  may  result  from  dislocation  glide,  deformation  twinning, 
and/or  sliding  of  micro-cracks.  For  uniaxial  strain  compression,  mode  I  crack  opening  is  omitted,  and  inelasticity  is  isochoric 
(sg  ■  mg  =  0  =>•  /  =  1 ).  From  geometry  of  the  problem,  all  n  systems  experience  the  same  resolved  shear  stress  x  =  Ta.  Nor¬ 
malizing  (3.14)  by  shear  modulus  (2.82)  or  (2.83),  the  prescribed  yield  criterion  in  the  plastically  deforming  regime  is 

T/Go=g(£)/Go=go-X£;  £  =  ny.  (3.22) 


Here,  g0  is  dimensionless  initial  shear  strength  at  the  HEL,  dependence  of  strength  g  =  ga  on  temperature  is  omitted,  and 
internal  state  variable  £  is  the  total  slip  on  all  n  systems.  Material  constant  X  >  0  accounts  for  loss  of  shear  resistance  as 
fracture  ensues  in  conjunction  with  dislocation  slip  or  twinning.  For  the  brittle  solids  considered  here,  such  strength  loss 
far  exceeds  any  hardening  that  might  arise  from  dislocation  accumulation.  In  U  =  U  of  (3.15),  it  is  assumed  k  =  0  since  prior 
analysis  (Clayton,  2009,  2010a,  2010c,  2011b)  confirms  that  stored  energy  of  lattice  defects  should  be  negligible  relative  to 
plastic  dissipation  in  strong  ceramics  under  uniaxial  compression,  insignificantly  affecting  the  predicted  thermomechanical 
response. 

Crystalline  quartz  is  known  to  have  a  large  number  of  possible  cleavage  planes  with  very  similar  surface  energy,  and  is 
also  prone  to  conchoidal  fracture  (Schultz,  Jensen,  &  Bradt,  1994),  i.e.,  curved  failure  surfaces  not  confined  to  any  particular 
plane,  as  observed  in  glass.  Dislocation  slip  (Clayton,  Chung,  Grinfeld,  &  Nothwang,  2008)  and  Dauphne  twinning  (Barton  & 
Wenk,  2009)  which  may  occur  at  high  temperatures  are  omitted  in  the  present  description  of  quartz  under  shock  compres¬ 
sion.  Therefore,  inelastic  deformation  consisting  of  modell/HI  cracks  on  numerous,  possibly  curved  surfaces  in  quartz  is  mod¬ 
eled  as  isochoric  and  isotropic: 


Fp(y) 


'  1  -  y  0  0 

o  (i-yr1/2  o 

0  0  (1  -yr1/2 


(3.23) 


Yield  criterion  (3.22)  applies  with  n  =  1  and  maximum  shear  stress  x  =  Jja^  -  cr3 1  used  in  place  of  t“,  where  (oq,  a2,  oq)  are 
principal  Cauchy  stress  components. 

Assume  that  HEL  shock  stress  PH  is  known  from  experimental  data.  Then  the  upstream  (HEL)  state  is  fully  determined  by 
the  solution  in  Section  2.2,  specifically  (2.73)-(2.76)  and  (2.77).  Specifically,  e  =  Iny  =  InF  is  decreased  incrementally  until 
(2.74)  reaches  PH,  at  which  point  F  =F+  and  U  =  U+.  Given  total  deformation  F~  and  slip  variable  y,  thermoelastic  deforma¬ 
tion  in  behind  the  plastic  shock  is  known  from  FE  =  F(F~)Fp~1(y).  Internal  energy,  axial  shock  stress,  and  shear  stress  can 
then  be  written  in  the  form 


U~  =  U~(F~, y,9]~),  P- =p-(F-,y,>r),  x  =  x(F~,y,  ij~).  (3.24) 

Let  F~  =  (Y-)~  be  prescribed  as  the  load  parameter.  Then  energy  balance  (3.19)  and  yield  criterion  (3.22)  comprise  two  cou¬ 
pled  algebraic  equations  that  can  be  solved  simultaneously  for  y  and  i/~: 

U-(F-,y,tj-)  ~  U+  =  i[P-(F-,y,fr)  +Ph][F+  -F-], 

?(F-,y,>r)/G0  =  g0-ny- 


(3.25) 

(3.26) 
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To  obtain  Hugoniot  stress  versus  volume  curves  reported  later  in  Section  3.3,  (3.26)  and  (3.25)  are  solved  simultaneously  for 
y  and  i/-  as  F~  is  decreased  incrementally  below  ^  from  the  HEL  state.  With  shock  stress  computed  from  the  second  of  (3.24), 
plastic  shock  velocity  25  and  downstream  particle  velocity  tr  can  be  obtained  from  the  Hugoniot  equations  for  mass  and 
momentum  conservation,  (3.16)  and  (3.17),  leading  to  (Molinari  &  Ravichandran,  2004) 

35  =  {(P--PH)/[P0(F+-F-)]}1/2,  !r  =  u+-:D(F--F+).  (3.27) 

The  downstream  state  is  now  fully  known.  The  above  procedure  signifies  the  first  known  (reported)  analytical  solution  for 
shock  compression  of  anisotropic  nonlinear  thermoelastic-plastic  single  crystals,  albeit  restricted  to  planar  shocks  in  highly 
symmetric  crystal  orientations.  Previous  analytical  solutions  have  been  restricted  to  linear  isentropic  elasticity  (Johnson, 
1972;  Johnson,  1974;  Johnson  et  al.,  1970)  in  anisotropic  crystals  or  to  finite  strain  isotropic  elastoplasticity  (Perrin  & 
Delannoy-Coutris,  1983). 

3.3.  Results:  elastic-plastic  solution 

Elastic-plastic  solutions  of  Section  3.2  are  now  applied  to  three  single  crystals:  Z-cut  sapphire,  X-cut  diamond,  and  X-  and 
Z-cut  quartz.  Inelastic  deformation  mechanisms  contributing  to  cumulative  shear  y  and  inelastic  properties  n.g0,  and  % 
(respectively  the  number  of  shear  systems,  initial  dimensionless  shear  strength,  and  rate  of  strength  reduction  with  cumu¬ 
lative  inelastic  shear)  are  listed  in  Table  3.  Also  shown  is  the  predicted  value  of  cumulative  shear  yr  at  which  the  Hugoniot 
stress  approaches  the  hydrostatic  pressure-volume  curve  of  (A.2)  given  in  the  Appendix  A  (P  — >  p  as  y  — >  yF).  Particular  en¬ 
tries  of  Table  3  are  discussed  in  detail  for  each  material  in  conjunction  with  predicted  Hugoniot  stress-volume  curves  in  what 
follows  next. 

For  Z-cut  sapphire,  preferred  deformation  mechanisms  are  rhombohedral  (R-plane)  twinning  and  R-plane  fracture,  which 
occur  equally  on  three  {101  2} (1011)  shear  systems  (n  =  3).  Previous  analysis  and  experiments  have  confirmed  that 
R-plane  twinning  is  preferred  over  all  other  glide  and  twinning  mechanisms  for  this  orientation  (Clayton,  2009;  Fuller, 
Winey,  &  Gupta,  2013).  Fracture  is  also  most  likely  to  occur  on  these  planes,  which  have  relatively  low  surface  energy 
(Schultz  et  al.,  1994);  dynamic  R-plane  cleavage  has  been  confirmed  by  high  speed  photography  in  impact  experiments 
(McCauley,  Strassburger,  Patel,  Paliwal,  &  Ramesh,  2013).  Initial  yield  strength  for  twinning  at  the  HEL  is  predicted  as 
g  «  0.044Go  =  7.3  GPa,  near  the  upper  bound  of  ranges  1  <g  <  8  GPa  quoted  for  R-twinning  from  previous  analyses  of 
indentation  (Tymiak  &  Gerberich,  2007)  and  uniaxial  strain  (Clayton,  2009).  Noting  that  the  characteristic  twinning  shear 
for  this  twin  system  is  yT  =  0.202  (Clayton,  2009),  the  maximum  volume  fraction  of  twinned  material  for  each  system  is 
«  0.3.  Predicted  Hugoniot  stress  is  compared  with  experimental  data  (Graham  &  Brooks,  1971)  in  Fig.  4.  The  HEL 


Table  3 

Inelastic  deformation  mechanisms  and  properties  for  single  crystals. 


Material 

Orientation 

Mechanism(s) 

n 

go 

X 

7f 

Sapphire 

Z 

{101 2}<1 0  11)  twinning  and  cleavage 

3 

0.044 

0.25 

0.06 

Diamond 

X 

(1 1 1}(1 1 0)  glide  and  cleavage 

8 

0.041 

0.10 

0.04 

X 

{111}  (2 11)  glide  and  cleavage 

4 

0.047 

0.10 

0.07 

Quartz 

X 

conchoidal  fracture 

1 

0.047 

0.33 

0.10 

Z 

conchoidal  fracture 

1 

0.087 

0.50 

- 

VI V0 


Fig.  4.  Nonlinear  elastic-plastic  solution  and  experimental  data  (Graham  &  Brooks,  1971)  for  Z-cut  sapphire:  shock  stress  vs.  volume  ratio. 
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Fig.  5.  Nonlinear  elastic-plastic  solutions  and  experimental  data  (Lang  &  Gupta,  2010;  Kondo  &  Ahrens,  1983;  Pavlovskii,  1971)  for  diamond:  shock  stress 
vs.  volume  ratio. 


O 

q: 


viv0  v/v0 


(a) 


(b) 


Fig.  6.  Nonlinear  elastic-plastic  solutions  and  experimental  data  (Fowles,  1967)  for  quartz:  (a)  shock  stress  vs.  volume  ratio,  X-cut  (b)  shock  stress  vs. 
volume  ratio,  Z-cut. 


corresponds  to  ^  a:  0.96,  at  which  point  twinning  initiates  according  to  the  model,  in  the  plastic  regime,  Hugoniot  stress  col¬ 
lapses  to  the  hydrostat  of  (A.2)  for  ^  <  0.88.  Since  strength  decreases  >  0),  fracture  on  k-planes  accompanies  twinning. 

For  X-cut  diamond,  two  possible  preferred  deformation  mechanisms  are  considered  separately:  (i)  octahedral  slip  and 
octahedral  plane  fracture  which  occur  equally  on  eight  {1 1 1  }(1 1 0)  shear  systems  (n  =  8)  and  (ii)  slip  and  octahedral  plane 
fracture  which  occur  equally  on  four  {111}  (2 11)  shear  systems  (n  =  4).  These  two  sets  of  systems  have  been  suggested  else¬ 
where  (Lang  &  Gupta,  2010  and  references  therein)  for  diamond,  the  first  being  reported  more  often.  Fracture  is  most  likely 
to  occur  on  { 1 1 1 }  planes,  which  have  low  surface  energy  (Schultz  et  al.,  1 994)  and  low  cleavage  strength  predicted  from  first 
principles  (Telling,  Pickard,  Payne,  &  Field,  2000).  Initial  yield  strength  for  slip  at  the  HEL  is  predicted  as 
g  m  0.041  Go  =  22  GPa  for  (T 1  0)  systems  and  g  «  0.047G0  =  25  GPa  for  (211)  systems.  These  predictions  are  within  ranges 
19  <  g  <  30  GPa  for  (1 1  0)  and  23  <  g  <  35  GPa  for  (21 1)  estimated  in  Lang  and  Gupta  (2010)  using  isentropic  Lagrangian 
elasticity.  Using  isotropic  elasticity,  a  shear  strength  of  30  GPa  has  also  been  estimated  from  shock  data  (Kondo  &  Ahrens, 
1983).  Predicted  Hugoniot  stress  is  compared  with  experimental  data  (Kondo  &  Ahrens,  1983;  Lang  &  Gupta,  2010; 
Pavlovskii,  1971)  in  Fig.  5.  While  shock  experiments  were  along  [100]  directions  in  Lang  and  Gupta  (2010),  Kondo  and  Ahrens 
(1983),  and  Pavlovskii  (1971),  shock  experiments  in  Kondo  and  Ahrens  (1983)  were  not  along  pure  mode  directions,  so  mod¬ 
el  predictions  are  only  approximately  representative  of  the  latter  data.  In  the  model,  the  representative  HEL  corresponds  to 
4L  rs  0.94,  towards  the  upper  end  of  the  range  of  values  reported  in  Lang  and  Gupta  (2010)  and  consistent  with  experiments 
of  Kondo  and  Ahrens  (1983).  In  the  plastic  regime,  Hugoniot  stress  collapses  to  the  hydrostat  of  (A.2)  for  X.  <  0.82  for  (T  1  0) 
systems  or  ^  5  0.80  for  (211)  systems.  Strength  decreases  at  the  same  rate  (/  =  0.10)  for  each  model  prediction,  meaning 
fracture  on  octahedral  planes  accompanies  dislocation  glide  in  each  case.  While  both  sets  of  slip  systems  enable  accurate 
representation  of  experimental  data,  the  present  calculations  suggest  that  {1 1 1  }(1 10)  slip  should  be  more  likely  to  occur 
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(a) 


(b)  (c) 


Fig.  7.  Nonlinear  elastic-plastic  solutions:  (a)  cumulative  slip  or  inelastic  shear  (b)  temperature  (c)  entropy. 


than  {111}(211)  slip  since  g0  is  lower  in  the  former.  More  slip  accumulates  on  each  individual  system  in  the  latter  set  (i.e., 
larger  yF  in  Table  3  for  (2 1 1))  since  there  are  fewer  systems  (n  =  4)  for  this  set. 

For  X-  and  Z-cut  quartz,  conchoidal  fracture  or  simultaneous  cleavage  of  a  large  number  of  planes  (Schultz  et  al ,  1994)  is 
modeled  using  the  isotropic  representation  of  inelastic  deformation  in  (3.23).  Yield  is  permitted  to  be  anisotropic  through 
prescription  of  different  values  of  g0  and  '/  depending  on  orientation.  This  approach  enables  successful  fitting  to  Hugoniot 
data  (Fowles,  1967)  in  Fig.  6,  but  is  not  fully  predictive  since  experimental  data  are  insufficient  to  parameterize  shear 
strength  for  arbitrary  crystal  orientations.  Initial  yield  strength  for  fracture  at  the  HEL  is  computed  asg«  0.047G0  =  2.3  GPa 
forX-cut  quartz  andg  «  0.087G0  =  4.2  GPa  for  Z-cut.  In  the  model,  the  representative  HEL  corresponds  to  X-  ks  0.94  forX-cut 
and  41  Kb  0.92  for  Z-cut.  In  the  plastic  regime,  Hugoniot  stress  nears  the  hydrostat  of  (A.2)  for  X.  <  0.84  for  X-cut  quartz  in 
Fig.  6(a),  wherein  y  rs  yF  =  0.1.  In  Fig.  6(b),  data  from  Fowles  (1967)  remain  above  the  logarithmic  hydrostat  to  volumetric 
compression  in  excess  of  20%.  For  this  reason,  no  value  of  yF  is  listed  for  Z-cut  quartz  in  Table  3,  and  %  — >  0  for  y  >  0.092  is 
imposed  in  calculations  for  this  orientation. 

Theoretical  predictions  of  cumulative  slip,  temperature,  and  entropy  are  shown  for  all  crystal  types  in  respective  Fig.  7(a)-(c). 
All  three  variables  are  closely  related  since  entropy  and  temperature  increase  in  conjunction  with  energy  dissipated  by  inelas¬ 
tic  deformation  in  the  plastic  regime.  Temperature  rise  and  entropy  production  are  comparatively  small  in  the  elastic  regime 
[Fig.  7(b)  and  (c)]  as  noted  already  in  Section  2.3  and  Table  2.  Predicted  temperature  rise  and  normalized  entropy  production 
are  largest  for  diamond  because  it  has  by  far  the  largest  absolute  shear  strength  g  and  the  smallest  c0  (Table  1 )  of  the  three 
crystals. 

4.  Conclusions 

Finite  strain  thermoelastic  and  thermoelastic-plastic  constitutive  theories  incorporating  a  logarithmic  elastic  strain 
measure  referred  to  material  coordinates  have  been  developed  for  anisotropic  single  crystals  and  applied  towards  shock 
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Fig.  8.  Logarithmic  EOS  predictions  and  experimental  data  (Hart  &  Drickamer,  1965;  Occelli  et  al„  2003;  Kimizuka  et  al„  2007)  for  hydrostatic  compression 
of  sapphire,  diamond,  and  quartz. 


compression  problems.  New  analytical  solutions  have  been  derived  for  planar  elastic  and  plastic  shocks  along  directions  of 
high  symmetry.  Predictions  have  been  analyzed  for  sapphire,  diamond,  and  quartz  single  crystals.  In  the  elastic  regime,  supe¬ 
rior  accuracy  of  the  proposed  logarithmic  theory  over  existing  Lagrangian  and  Eulerian  theories  has  been  demonstrated  for 
sapphire,  diamond,  and  Z-cut  quartz.  In  the  plastic  regime,  constitutive  relations  and  necessary  parameters  have  been  devel¬ 
oped  and  determined  that  describe  experimental  Hugoniot  data,  including  the  loss  of  shear  strength  for  sapphire,  diamond, 
and  X-cut  quartz  observed  at  shock  stresses  exceeding  the  HEL  and  the  corresponding  collapse  of  the  Hugoniot  stress  to  the 
hydrostat.  Accuracy  of  the  latter  described  by  a  logarithmic  pressure-volume  EOS  has  been  demonstrated.  Precise  values  of 
rhombohedral  twinning  resistance  in  sapphire  and  octahedral  slip  resistance  in  diamond  have  been  calculated;  these 
compare  favorably  with  coarse  estimates  given  elsewhere.  Model  results  for  quartz  suggest  that  fracture  criteria  are  highly 
orientation  dependent  despite  the  availability  of  numerous  low-energy  planes  for  potential  cleavage  and  experimental  evi¬ 
dence  of  conchoidal  fracture. 


Appendix  A.  Hydrostatic  compression 


Considered  here  is  the  pressure-volume  equation-of-state  (EOS)  derived  from  the  logarithmic  nonlinear  elasticity  theory 
of  Sections  2.1  and  3.1.  Assume  adiabatic  conditions,  isochoric  plastic  deformation  (if  any)  such  that  J  =  detF  =  JE  consists 
only  of  thermoelastic  volume  change,  and  a  hydrostatic  stress  state  such  that  a  =  —pi,  with  p  the  Cauchy  pressure.  Then 
balance  of  energy  (2.7)  degenerates  to 

U  =Ja  :  Vv  =  -JpV  ■  v  =  -pj.  (A.l) 


For  materials  with  less  than  cubic  symmetry,  deviatoric  components  of  deformation  rate  may  be  nonzero,  but  only  volume 
changes  contribute  to  internal  energy.  Pressure  obeys,  for  logarithmic  nonlinear  elasticity,  the  following  EOS  (Poirier  & 
Tarantola,  1998): 


P  =  -(d0/dj)  |A„=0 


-Bo[(ln  ])/J] 


2)  In  J 


(A.2) 


where  B0  is  the  isentropic  bulk  modulus  and  B(,  the  pressure  derivative  of  the  bulk  modulus  in  the  unstressed  reference  state. 
For  isothermal  hydrostatic  loading,  free  energy  and  isothermal  constants  replace  internal  energy  and  isentropic  constants  in 
(A.2).  For  2  ^  Bj,  ^  5  and  §  ^  1,  differences  in  pressure  predicted  by  logarithmic  EOS  (A.2)  and  the  second-order  Birch- 

Murnaghan  (B-M)  EOS  associated  with  Eulerian  elasticity  (Clayton,  2013;  Jeanloz,  1989)  are  insignificant,  but  at  very  high 
compression  -  as  occurring  in  the  Earth’s  interior  for  example  -  the  logarithmic  EOS  becomes  more  realistic  than  the  B-M 
EOS  (Poirier  &  Tarantola,  1998). 

Predictions  of  (A.2)  are  compared  with  experimental  hydrostatic  compression  data  on  sapphire  (Hart  &  Drickamer,  1965), 
diamond  (Occelli,  Loubeyre,  &  LeToullec,  2003),  and  quartz  (Kimizuka,  Ogata,  Li,  &  Shibutani,  2007)  in  Fig.  8.  Pressures  are 
limited  to  those  under  which  phase  transformations  might  occur.  Close  agreement  with  all  data  is  evident  for  >  0.88.  Pre¬ 
dictions  are  made  using  isentropic  bulk  moduli  from  (2.80)  or  (2.81 )  and  pressure  derivatives  obtained  from  ultrasonic  data 
(Gieske  &  Barsch,  1968;  McSkimin  &  Andreatch,  1972;  McSkimin  et  al.,  1965),  all  listed  in  Table  1,  i.e.,  the  EOS  parameters 
were  obtained  independently  and  are  not  fit  directly  to  the  high  pressure  data  shown.  The  EOS  predictions  tend  to  exceed  the 
data  for  diamond  and  quartz  for  ^  <  0.88;  the  discrepancy  is  due  in  part  to  use  of  isentropic  elastic  constants  (more  appro¬ 
priate  for  comparison  with  shock  data  in  Section  3.3)  rather  than  isothermal  constants. 
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